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Abstract 

We propose a new algorithm for sampling the A-body density |^(R)| 2 / L 3JV |^| 2 in the Vari- 
ational Monte Carlo (VMC) framework. This algorithm is based upon a modified Ricci-Ciccotti 
discretization of the Langevin dynamics in the phase space (R, P) improved by a Metropolis accep- 
tation/rejection step. We show through some representative numerical examples (Lithium, Fluorine 
and Copper atoms, and phenol molecule), that this algorithm is superior to the standard sampling 
algorithm based on the biased random walk (importance sampling). 

PACS numbers: 
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I. INTRODUCTION 



Most quantities of interest in quantum physics and chemistry are expectation values of 
the form 

1 ' 

where O is the self-adjoint operator (the observable) associated with a physical quantity O 
and ^ a given wave function. 

For iV-body systems in the position representation, \I/ is a function of 3N real variables 
and 

1 [0^](R)^(R)*dR 



(jS_Ojr) 



|#(R)| 2 dR 



(2) 



High-dimensional integrals are very difficult to evaluate numerically by standard integration 
rules. For specific operators O and specific wave functions e.g. for electronic Hamiltonians 
and Slater determinants built from Gaussian atomic orbitals, the above integrals can be 
calculated analytically. In some other special cases, (j2J) can be rewritten in terms of integrals 
on lower-dimensional spaces (typically M 3 or M 6 ). 

In the general case however, the only possible way to evaluate (J2J) is to resort to stochastic 
techniques. The VMC method |jj consists in remarking that 



L (R) |^(R)| 2 dR 

(3) 



|^(R)| 2 dR 



with L (R) = [0tf](R)/tf(R), hence that 

N ' ' n=l 

where (R n ) n >i are points of M. 3N drawn from the probability distribution |vI/(R)| 2 / J R3N |\[>| 2 . 

The VMC algorithms described in the present article are generic, in the sense that they 
can be used to compute the expectation value of any observable, for any iV-body system. 
In the numerical example, we will however focus on the important case of the calculation of 
electronic energies of molecular systems. In this particular case, the expectation value to be 
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computed reads 

f ^WI^)l 2 dR 



(5) 



|^(R)| 2 dR 



where the scalar field £x(R) = [i^ l I / ](R)/ l I / (R) is called the /oca/ energy. Remark that if iff 
is an eigenfunction of associated with the eigenvalue E, -El(R) = E for all R. Most often, 
VMC calculations are performed with trial wave functions iff that are good approximations 
of some ground state wave function \lV Consequently, -El(R) usually is a function of low 
variance (with respect to the probability density |\1/(R)| 2 / J R3N \iff\ 2 )- This is the reason why, 
in practice, the approximation formula 

(iff Hiff) 1 

{ ' J ~-V£ L (R") (6) 
(iff, iff) L ^ K } 

1 ' n=l 

is fairly accurate, even for relatively small values of L (in practical applications on realistic 
molecular systems L ranges typically between 10 6 and 10 9 ). 

Of course, the quality of the above approximation formula depends on the way the points 



(R n ) n >i are generated. In section III Bl we describe the standard sampling method cur- 
rently used for VMC calculations. It consists in a biased (or importance sampled) random 
walk in the configuration space (also called position space) M. 3N corrected by a Metropolis 
acceptation/rejection procedure. In section Hi C\ we introduce a new sampling scheme in 
which the points (R n ) n >i are the projections on the configuration space of one realization of 
some Markov chain on the phase space (also called position-momentum space) M 3N x M. 3N . 
This Markov chain is obtained by a modified Langevin dynamics, corrected by a Metropolis 
acceptation/rejection procedure. 

Finally, some numerical results are presented in section ITTTl Various sampling algorithms 
are compared and it is demonstrated on a bench of representative examples that the algo- 
rithm proposed here based on the modified Langevin dynamics is the most efficient one (the 
mathematical criteria for measuring the efficiency will be made precise below) . 

Before turning to the technical details, let us briefly comment on the underlying moti- 
vations of our approach. The reason why we have introduced a (purely fictitious) Langevin 
dynamics in the VMC framework is twofold: 

• first, sampling methods based on Langevin dynamics turn out to outperform those 
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based on biased random walks in classical molecular dynamics (see [4] for a quantitative 
study on carbon chains); 

• second, a specific problem encountered in VMC calculations on fermionic systems is 
that the standard discretization of the biased random walk (Euler scheme) does not 
behave properly close to the nodal surface of the trial wave function This is due to 
the fact that the drift term blows up as the inverse of the distance to the nodal surface: 
if a random walker gets close to the nodal surface, the drift term repulses it far apart 
in a single time step. As demonstrated in 0,0], it is possible to partially circumvent 
this difficulty by resorting to more clever discretization schemes. Another strategy 
consists in replacing the biased random walk by a Langevin dynamics: the walkers 
then have a mass (hence some inertia) and the singular drift does not directly act on 
the position variables (as it is the case for the biased random walk), but indirectly via 
the momentum variables. The undesirable effects of the singularities are thus expected 
to be damped down. 

II. DESCRIPTION OF THE ALGORITHMS 
A. Metropolis algorithm 

The Metropolis algorithm is a general purpose sampling method, which combines the 
simulation of a Markov chain with an acceptation/rejection procedure. 

In the present article, the underlying state space is either the configuration space R 3N or 
the phase space R 3N x R 3N = R 6N . Recall that a Markov chain on R d is characterized by 
its transition kernel p. It is by definition the non-negative function of R d x B(R d ) (B(R d ) is 
the set of all the Borel sets of R d ) such that, if X e R d and B 6 B{R d ), the probability for 
the Markov chain to lay in B at step n + 1 if it is at X at step n is p(X, B). The transition 
kernel has a density with respect to the Lebesgue measure if for any X 6 R d , there exists a 
non-negative function fx € L 1 (R d ) such that 



The non-negative number /x(X') is often denoted by T(X — > X') and the function T 




(7) 



R d x R 



.+ is called the transition density. 
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Given a Markov chain on M d with transition density T and a positive function / G L 1 { 
the Metropolis algorithm consists in generating a sequence (X™) ne N of points in R d starting 
from some point X° G M. d according to the following iterative procedure: 

• propose a move from X" to X n+1 according to the transition density T(X n — > X n+1 ); 

• compute the acceptance rate 

A(X» X»«) = min ( /(X" + ')r(X- + - r X , 

\ /(X")T(X n X" +1 ) 

• draw a random variable £7 n uniformly distributed in [0, 1] ; 

- if U n < A(X n -> X n+1 ), accept the move: X n+1 = X n+1 ; 

- if U n > A(X n -> X n+1 ), reject the move: X n+1 = X n . 



It is not difficult to show (see p| for instance) that, for a very large class of transition densities 
T, the points X n generated by the Metropolis algorithm are asymptotically distributed 
according to the probability density /(X)/ j Rd f. On the other hand, the practical efficiency 
of the algorithm crucially depends on the choice of the transition density (i.e. of the Markov 
chain) . 



B. Random walks in the configuration space 

In this section, the state space is the configuration space M 3N and / = |\1/| 2 , so that the 
Metropolis algorithm actually samples the probability density |\I/(R)| 2 / f R3N |\&| 2 . 



1. Simple random walk 

In the original paper of Metropolis et al, the Markov chain is a simple random walk: 

R n+1 = R" + AR U n (8) 

where AR is the step size and U n are independent and identically distributed (i.i.d.) random 
vectors drawn uniformly in the 3iV-dimensional cube K = [— 1,1} 3N . The corresponding 
transition density is T(R — > R') = 2~ 3N xk ((R — ~R')/AR) where xk is the characteristic 
function of the cube K (note that in this particular case, T(R — > R') = T(R' — > R)). 
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2. Biased random walk 



The simple random walk is far from being the optimal choice: it induces a high rejection 
rate, hence a large variance. A variance reduction technique usually referred to as the 
importance sampling method, consists in considering the so-called biased random walk or 



over-damped Langevin dynamics [11]: 

dR(t) = V[log |tf |](R(t))dt + dW(t), (9) 

where W(£) is a 3iV-dimensional Wiener process. Note that |\I/| 2 is an invariant measure of 
the Markov process (JHJ), and, better, that the dynamics (jSJ) is in fact ergodic and satisfies a 
detailed balance property The qualifier ergodic means that for any compactly supported 
continuous function g : R 3N — > R, 

s(R)|vf(R)| 2 dR 



lim If g(R(t))dt = ^- . (10) 

T ^ +ooTJo / |M>(R)| 2 dR 

The detailed balance property reads 

|M>(R)| 2 T At (R^R') = |M/(R')| 2 T A4 (R'^R) (11) 

for any At > 0, where T At (R — > R') is the probability density that the Markov process Q 
is at R' at time t + At if it is at R at time t. These above results are classical for regular, 

n, 

positive functions ty, and have been recently proven for fermionic wave functions [lj] (in the 
latter case, the dynamics is ergodic in each nodal pocket of the wave function 

Note that if one uses the Markov chain of density T A t(R — > R') in the Metropolis algo- 
rithm, the acceptation/rejection step is useless, since due to the detailed balance property, 
the acceptance rate always equals one. 

The exact value of T A <(R — > R') being not known, a discretization of equation Q with 
Euler scheme, is generally used 

R n+1 = R n + AtV[log|^|](R") + AW n (12) 

where AW n are i.i.d. Gaussian random vectors with zero mean and covariance matrix 
At I 3 n {Iw is the identity matrix). The Euler scheme leads to the approximated transition 
density 
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(R ^ R) -(2^)W^ eXP l v 2Ai J (13) 

The time discretization introduces the so-called time-step error, whose consequence is that 
(fT2|) samples |\I/(R)| 2 / f R3N |\l/| 2 | only approximately. Note that the Metropolis accepta- 
tion/rejection procedure perfectly corrects the time-step error. In the limit At — ► 0, the 
time-step error vanishes and the acceptation/rejection procedure is useless. 

This sampling method is much more efficient than the Metropolis algorithm based on the 
simple random walk, since the Markov chain (fT2|) does a large part of the work (it samples 
a short time-step approximation of ^(R)) 2 / f \^ 2 \), which is clearly not the case for the 
simple random walk. 

The standard method in VMC computations currently is the Metropolis algorithm based 
on the Markov chain defined by l|12J) . For refinements of this method, we refer to 

C. Random walks in the phase space 

In this section, the state space is the phase space M. 3N x M. 3N . Let us emphasize that the 
introduction of momentum variables in nothing but a numerical artifice. The phase space 
trajectories that will be dealt with in this section do not have any physical meaning. 

1. Langevin dynamics 

The Langevin dynamics of a system of N particles of mass m evolving in an external 
potential V reads 

' dR(t) = ±P(t)dt, 
dP(t) = -W(R(t))dt-jP(t)dt + adW(t). 

As above, R(i) is a 3iV-dimensional vector collecting the positions at time t of the iV par- 
ticles. The components of the 3iV-dimensional vector P(£) are the corresponding momenta 
and W(t) is a 3iV-dimensional Wiener process. The Langevin dynamics can be considered 
as a perturbation of the Newton dynamics (for which 7 = and a — 0). The magnitudes a 
and 7 of the random forces crdW(i) and of the drag term — 7P(t)di are related through the 
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fluctuation-dissipation formula 



°>= 2 -f, (IB) 

where (3 is the reciprocal temperature of the system. Let us underline that in the present 
setting, (3 is a numerical parameter that is by no means related to the physical temperature 
of the system. It can be checked (at least for regular potentials V) that the canonical 
distribution 

dII(R, P) = Z-V^^dRdP (16) 
is an invariant probability measure for the system, Z being a normalization constant and 



IP 
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£T(P,R) = V(R) + ^ (17) 

being the Hamiltonian of the underlying Newton dynamics. In addition, the Langevin 
dynamics is ergodic (under some assumptions on V). Thus, choosing 

(3 = 1 and V = -log|^| 2 , (18) 

the projection on the position space of the Langevin dynamics samples |\I/(R)| 2 / J |\I/| 2 . On 
the other hand, the Langevin dynamics does not satisfy the detailed balance property. We 
will come back to this important point in the forthcoming section. 

In this context, the parameters m and 7 (ex being then obtained through (|TK|) ) should be 
seen as numerical parameters to be optimized to get the best sampling. We now describe 
how to discretize and apply a Metropolis algorithm to the Langevin dynamics |JTJJ, in the 
context of VMC. 



2. Time discretization of the Langevin dynamics 

Many discretization schemes exist for Langevin dynamics. In order to choose which 
algorithm is best for VMC, we have tested four different schemes available in the litera- 
ture 0HHQ, with p.— , = 1, 7 = 1 and m = 1. 0. be „ chmarli SyS - 
tem is a Lithium atom, and a single determinantal wave function built upon Slater-type 
atomic orbitals, multiplied by a Jastrow factor. We turn off the acceptation/rejection step 
in these preliminary tests, since our purpose is to compare the time-step errors for the var- 
ious algorithms. From the results displayed in table HJ one can see that the Ricci-Ciccotti 
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algorithm [LJ] is the method which generates the smallest time-step error. This algorithm 
reads 



where G n are i.i.d. Gaussian random vectors with zero mean and variance cr 2 / 3Ar with 
a 2 = 2f±At. 

It can be seen from Table I that the Ricci-Ciccotti algorithm also outperforms the biased 
random walk (JTJJ, as far as sampling issues are concerned. In the following, we shall therefore 
use the Ricci-Ciccotti algorithm. 

3. Metropolized Langevin dynamics 

The discretized Langevin dynamics does not exactly sample the target distribution IT, but 
rather from some approximation IIa* of IT. It is therefore tempting to introduce a Metropolis 
acceptation/rejection step to further improve the quality of the sampling. Unfortunately, 
this idea cannot be straightforwardly implemented for two reasons: 

• first, this is not technically feasible, for the Markov chain defined by (fT9| does not have 
a transition density. Indeed, as the same Gaussian random vectors G n are used to 
update both the positions and the momenta, the measure p((R™, P n ), ■) is supported 
on a 3./V-dimensional submanifold of the phase space M? N x IR 3Ar ; 

• second, leaving apart the above mentioned technical difficulty, which is specific to the 
Ricci-Ciccotti scheme, the Langevin dynamics is a priori not an efficient Markov chain 
for the Metropolis algorithm, for it does not satisfy the detailed balance property. 

Let us now explain how to tackle these two issues, starting with the first one. 

To make it compatible with the Metropolis framework, one needs to slightly modify the 
Ricci-Ciccotti algorithm. Following 0, QiJ, we thus introduce i.i.d. correlated Gaussian 




>Ai/2 + m [ - W(R n ) At + G n ] e"^ 4 , 
[W(R n ) + Vy(R" +1 )]e-^/ 2 + G n e^ At / 2 , 



(19) 
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vectors (G^, G^) (1 < % < 3N) such that: 

At 



Setting G™ 
reads 



(3 my 
m 



3 _ 4e -7At + e -2 7 Ai 



7 At 



Cl2 



" ' - and Go 



0102 
(^l,t)l<»<3JV 



(20a) 
(20b) 

— ^. (20c) 

P7aicr 2 

(^2i)i<<<3Jv 3 the modified Ricci-Ciccotti algorithm 



3 - 7 At\2 



(21) 



R n+1 = R n + Aipn e - 7 At/2 _ |£ (< R n) e - 7 At/4 + Q„ 
pn+1 = pn e - 7 Ai _ At [ V y( R n) + VV(R n+1 )] e"^ 2 + G£. 

The above scheme is a consistent discretization of and the corresponding Markov chain 
does have a transition density, which reads (see Appendix) 



T^ RC ((R r \P n ) -> (R n+1 ,P n+1 

z-< 



d! = R 



cxp 



n+1 




+ 



di d 2 
-2c,o — ■ — 



R n - At- 



m 



. e -7At/2 + 



At 2 
2m 



W(R n )e~ 7A ' /4 , 



d 2 = P n+1 - p n e~^ At + -At [W(R n ) + W(R n+1 )]e" 7A ' /2 . 



-i 1,11 / \ / -• i - i (22a) 

2 (! - c i2) \ V °"i / V a 2 J <?i cr 2 J 

with 

(22b) 
(22c) 

Unfortunately, inserting directly the transition density (|22|) in the Metropolis algorithm leads 
to a high rejection rate. Indeed, if (R n , P") and (R n+1 ,P n+1 ) are related through (|2Tj) . 
T^ RC ((R n ,P n ) -> (R n+1 ,P n+1 )) usually is much greater than T^ RC ((R n+1 , P n+1 ) -> 
(R n ,P n )), since the probability that the random forces are strong enough to make the 
particle go back in one step from where it comes, is very low in general. This is related to 
the fact that the Langevin dynamics does not satisfy the detailed balance relation. 

It is however possible to further modify the overall algorithm by ensuring some micro- 
scopic reversibility, in order to finally obtain low rejection rates. For this purpose, we 
introduce momentum reversions. Denoting by y^" ngevm the transition density of the Markov 
chain obtained by integrating ([TH) exactly on the time interval [t,t + At], it is indeed not 
difficult to check (under convenient assumptions on V = — log | \P | 2 , that the Markov chain 
defined by the transition density 
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T Langevin ((R> p) ^ = ^Langevin ((R p) ^ (R? _ p/)) 

is ergodic with respect to II and satisfies the detailed balance property 

n(R, P) f^ ngevin ((R ? P ) _> (R, p')) = n(R', P') f^ ngcvin ((R', P') -> (R, P)) . (24) 

Replacing the exact transition density j^angcvm ^ approximation T^ RC , we now consider 
the transition density 

((R, P) - (R', P')^ RC ) = ^ RC ((R, P) - (R', -P')) • (25) 
The new sampling algorithm that we propose can be stated as follows: 

• Propose a move from (R n , P n ) to (R n+1 , P n+1 ) using the transition density T^ RC . In 
other words, perform one step of the modified Ricci-Ciccotti algorithm 

pn+l = R n + Atpn - 7 At/2 _ Aj^yy/pn) - 7 At/4 + Qn 

* m 2m \ / J- ' (26) 

pn+i = P n e - 7 At _ m [W(R n ) + W(R n+1 )] e -^ A */ 2 + G5. 
and set (R n+1 ,P n+1 ) = (R™ +1 , -P™ +1 ) 

• Compute the acceptance rate 



A((R n ,P n ) -> (R n+1 ,P n+1 )) 

/ n(R n+1 , P" +1 ) T^ t RC ((R n+1 , P" +1 ) -> (R n , P»)) \ 

- mm ^ ^ f^ RC ((Rn pn) _^ (Rn+1 > pn+1)) ' J ' 

• Draw a random variable £/ n uniformly distributed in (0, 1) and 

- if U n < A((R n ,P n ) -> (R n+1 ,P n+1 )), accept the proposal: (R n+1 ,P n+1 ) = 

(Rn+l > pn+l) j 

- if U n > A((R n , P n ) -> (R n+1 , P n+1 )), reject the proposal, and set (R" + \ P~" +1 ) = 
(R",P n ). 

• Reverse the momenta 

(R n+l ; pn+1) = (p^ _p" +1 ) (27) 
11 



Note that a momentum reversion is systematically performed just after the Metropolis 
step. As the invariant measure IT is left unchanged by this operation, the global algo- 
rithm (Metropolis step based on the transition density T^ RC plus momentum reversion) 
actually samples EL The role of the final momentum reversion is to preserve the underly- 
ing Langevin dynamics: while the proposals are accepted, the above algorithm generates 
Langevin trajectories, that are known to efficiently sample an approximation of the target 
density EI. Numerical tests seem to show that, in addition, the momentum reversion also 
plays a role when the proposal is rejected: it seems to increase the acceptance rate of the 
next step, preventing the walkers from being trapped in the vicinity of the nodal surface 
tf-^O). 

As the points (R n , P n ) of the phase space generated by the above algorithm form a sampling 
of II, the positions (R n ) sample (^(R)) 2 / f R3N |\I/| 2 and can therefore be used for VMC 
calculations. 



III. NUMERICAL EXPERIMENTS AND APPLICATIONS 



A. Measuring the efficiency 

A major drawback of samplers based on Markov processes is that they generate se- 
quentially correlated data. The effective number of independent observations is in fact 
L e ff = L/N C0TT , where N covr is the correlation length, namely the number of successive corre- 
lated moves. 

In the following applications, we provide estimators for the correlation length N C0TT and 
for the so-called inefficiency r] (see below), which are relevant indicators of the quality of the 
sampling. In this section, following Stedman et al. [l2|, we describe the way these quantities 
are defined and computed. 

The sequence of samples is split into Nb blocks of Lb steps, where the number Lb is 
chosen such that it is a few orders of magnitude higher than N COIZ . The empirical mean of 
the local energy reads 

, n b l b 

(El)^ Lb = wy - £ MR*). (28) 



i=i 
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The empirical variance over all the individual steps is given by 

N B L B 



1° Nb ' Lb ? = ttV £ (^( Ri ) - <^W) 2 (29) 

1=1 

and the empirical variance over the blocks by 

^I^fEK"^) , (30) 

where -Eb^ is the average energy over block z: 



iLf 



j={i-l)L B +l 

Following [l^l, we define the correlation length as 



E B,i = j- E ^( Ri )- ( 31 ) 



N cmr = Jim _ r hm _ L B \ B J (32) 



and the inefficiency r/ of the run as: 

V = lim lim L B [a^ fl ' ifl ] 2 . (33) 

Nb— *oo Lb^oo 

On the numerical examples presented below, the relative fluctuations of the quantities 

r ^B ' ^ B 1 2 

Lg |^ B ^ B j 2 and L s [cr^ B,is ] 2 become small for L B > 50 and iV B > 50. 

The definition of these two quantities can be understood as follows. Since L B 3> N C0TT and 
only Lb/N cott are independent samples among the samples in the block, the central limit 
theorem yields Esi — (-ExWp H — -, = where G l are i.i.d. normal random variables. 

y/ L B /N COII 

Thus, in the limit Nb — > oo and L# — > oo, we obtain that a B = Lg J N — which yields (|3~2l . 
The inefficiency rj is thus equal to N C0TT a 2 and is large if the variance is large, or if the 
number of correlated steps is large. 

Using this measure of efficiency, we can now compare the sampling algorithms (the simple 
random walk, the biased random walk and the Langevin algorithm) for various systems. In 
any case, a Metropolis acceptation/rejection step is used. We found empirically from several 
tests that convenient values for the parameters of the Langevin algorithm are 7=1 and 
m = Z 3 / 2 where Z is the highest nuclear charge among all the nuclei. For each algorithm, 
we compare the efficiency for various values of the step length, namely the increment AR 
in the case of the simple random walk, and the time-step At for the other two schemes. For 
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a given algorithm, simple arguments corroborated by numerical tests show that there exists 
an optimal value of this increment: for smaller (resp. for larger increments), the correlation 
between two successive positions increases since the displacement of the particle is small 
(resp. since many moves are rejected), and this increases the number of correlated steps 

-^corr • 

One can notice on the results (see tables IT U IIH [ IIV | lv|) that a large error bar corresponds 
to large values for N C0TT and rj. The quantities N C0TT and rj are a way to refine the measure 
of efficiency, since the same length of error bar may be obtained for different values of the 
numerical parameters. 

Let us now present some numerical tests. We compare the algorithms and parameters 
at a fixed computational cost. The reference values are obtained by ten times longer VMC 
simulations. The error bars given in parenthesis are 60% confidence intervals. We also 
provide the acceptance rate (denoted by A in the tables) and, when it is relevant, the mean 
of the length of the increment R™ +1 — R n over one time-step (denoted by (|AR|) in the 
tables) for the biased random walk and the Langevin dynamics. 

B. Atoms 

Lithium. The Lithium atom was chosen as a first simple example. The wave function 
is the same as for the benchmark system used for the comparison of the various Langevin 
schemes, namely a single Slater determinant of Slater-type basis functions improved by a 
Jastrow factor to take account of the electron correlation. The reference energy associated 
with this wave function is —7.47198(4) a.u., and the comparison of the algorithms is given in 
table HU The runs were made of 100 random walks composed of 50 blocks of 1000 steps. For 
the simple random walk, the lowest values of the correlation length and of the inefficiency 
are respectively 11.4 and 1.40. The biased random walk is much more efficient, since the 
optimal correlation length and inefficiency are more than twice smaller, i.e. 4.74 and 0.55. 
The proposed algorithm is even more efficient: the optimal correlation length is 3.75 and 
the optimal inefficiency is 0.44. 

Fluorine. The Fluorine atom was chosen for its relatively "high" nuclear charge (Z — 9) , 
leading to a timescale separation of the core and valence electrons. The wave function is a 
Slater-determinant with Gaussian-type basis functions where the Is orbital was substituted 
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by a Slater-type orbital, with a reference energy of —99.397(2) a.u. The runs were made of 
100 random walks composed of 100 blocks of 100 steps. The results are given in table lHIl For 
the simple random walk, the lowest values of the correlation length and of the inefficiency 
are respectively 15.6 and 282. The biased random walk is again twice more efficient than the 
simple random walk, for which the optimal correlation length and inefficiency are 7.4 and 
137. The Langevin algorithm is more efficient than the biased random walk: the optimal 
correlation length is 5.3 and the optimal inefficiency is 102. 

Copper. We can go even further in the timescale separation and take the Copper atom 
[Z = 29) as an example. The wave function is a Slater determinant with a basis of Slater- 
type atomic orbitals, improved by a Jastrow factor to take account of the electron correlation. 
The reference energy is —1639.2539(24). The runs were made of 40 random walks composed 
of 500 blocks of 500 steps. From table IIV[ one can remark that the Langevin algorithm is 
again more efficient than the biased random walk, since the optimal correlation length and 
inefficiency are respectively 28.7 and 4027, whereas using the biased random walk, these 
values are 51.0 and 5953. 

C. The phenol molecule 

The Phenol molecule was chosen to test the proposed algorithm because it contains three 
different types of atoms (H, C and O). The wave function here is a single Slater determinant 
with Gaussian-type basis functions. The core molecular orbitals of the Oxygen and Carbon 
atoms were substituted by the corresponding atomic Is orbitals. The comparison of the 
biased random walk with the Langevin algorithm is given in tabled The optimal correlation 
length using the biased random walk is 10.17, whereas it is 8.23 with our Langevin algorithm. 
The optimal inefficiency is again lower with the Langevin algorithm (544) than with the 
biased random walk (653). 

D. Discussion of the results 

We observe that on our numerical tests, the Langevin dynamics is always more efficient 
than the biased random walk. Indeed, we notice that: 

• The error bar (or iV C0IT , or rf) obtained with the Langevin dynamics for an optimal 
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set of numerical parameters is always smaller than the error bar obtained with other 
algorithms (for which we also optimize the numerical parameters). 

• The size of the error bar does not seem to be as sensitive to the choice of the numerical 
parameters as for other methods. In particular, we observe on our numerical tests that 
the value At = 0.2 seems to be convenient to obtain good results with the Langevin 
dynamics, whatever the atom or molecule. 
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APPENDIX A: DERIVATION OF THE TRANSITION PROBABILITY (|22|) 



The random vector (di,d 2 ) (defined by (I22bjl - (|22cjl ) is a Gaussian random vector and 
therefore admits a density with respect to the Lebesgue measure in M 67V . If, for 1 < i < 3N, 
we denote by d\^ (resp. d 2)i ) the components of d x (resp. d 2 ), we observe that the Gaussian 
random vectors {di >i) d 2i ) are i.i.d. Therefore, the transition probability T((R n ,P") — > 
(R n+1 ,P n+1 )) reads 



T((R n ,P n ) 



:r 



n+l TDra+l\ 



))=Z~ l ( P (d lji ,do i )) 



3A' 



(Al) 



where Z is a normalization constant and p denotes the density (in IR 2 ) of the Gaussian 
random vectors d 2} i). 

From equations (f2TJ) . one can see that 



pn A -2 

m 2m 
d 2)l = P? +1 - P" e - 7A * + l-At[ViV(K n ) + ViV(K n+1 )]e^ At / 2 , 



is a Gaussian random vector with covariance matrix T 



. Thus 



p(di,d 2 



^TrVdetr) 'exp ^-^(rf 1 ,d 2 )r- 1 (rf 1 ,d 2 ) T ^) 



(A2) 
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Since V 1 



J_ £12 

C12 1 



{ZED is easily obtained from l|XT]) ~ lfl2j) . 
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TABLE CAPTIONS 



• Table HJ Comparison of different discretization schemes for Langevin dynamics. The 
reference energy is -7.47198(4) a.u. 

• Table HU : The Lithium atom: Comparison of the Simple random walk, the Biased 
random walk and the proposed Langevin algorithm. The runs were carried out with 
100 walkers, each realizing 50 blocks of 1000 steps. The reference energy is -7.47198(4) 
a.u. 

• Table Mil : The Fluorine atom : Comparison of the Simple random walk, the Biased 
random walk and the proposed Langevin algorithm. The runs were carried out with 
100 walkers, each realizing 100 blocks of 100 steps. The reference energy is -99.397(2) 
a.u. 

• Table HXl : The Copper atom: Comparison of the Biased random walk with the pro- 
posed Langevin algorithm. The runs were carried out with 40 walkers, each realizing 
500 blocks of 500 steps. The reference energy is -1639.2539(24) a.u. 

• Table IVl: The Phenol molecule : Comparison of the Biased random walk with the pro- 
posed Langevin algorithm. The runs were carried out with 100 walkers, each realizing 
100 blocks of 100 steps. The reference energy is -305.647(2) a.u. 

TABLES 
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At BRW BBK [15J Force interpolation Splitting [1JJ Ricci & Ciccotti [13J 

0.05 -7.3758(316) -7.4395(246) -7.4386(188) -7.4467(137) -7.4576(07) 

0.005 -7.4644(069) -7.4698(015) -7.4723(015) -7.4723(015) -7.4701(20) 

0.001 -7.4740(007) -7.4728(013) -7.4708(017) -7.4708(017) -7.4696(17) 

0.0005 -7.4732(010) -7.4700(023) -7.4709(022) -7.4708(022) -7.4755(26) 

TABLE I: Comparison of Langevin algorithms and the biased random walk (BRW) 
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TABLE II: The Lithium atom 
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TABLE IV: The Copper atom 
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TABLE V: The Phenol molecule 
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